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We consider the design criteria of a chemical mixing device based on a chaotic flow, with an 
emphasis on the steady-state devices. The merit of a reactor, defined as the Q-factor, is related to 
the physical dimension of the device and the molecular diffusivity of the reactants through the local 
Lyapunov exponents of the flow. The local Lyapunov exponent can be calculated for any given flow 
field and it can also be measured in experimental situations. Easy-to-compute formulae are provided 
to estimate the Q-factor given either the exact spatial dependence of the local Lyapunov exponent 
or its probability distribution function. The requirements for optimization are made precise in the 
context of local Lyapunov exponents. 
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I. INTRODUCTION 

A major concern in chemical engineering is to react two or more species to form a new chemical product [Q. To 
increase production, one could boost the reaction rate by a catalyst which is an optimization on a microscopic level. 
One could also carry the reactants by a flow and maximize the diffusive effect that brings the chemicals together 
to react. Traditionally the flow is designed to be turbulent to achieve higher mixing effect. Started in the 1980s, 
non-turbulent but chaotic flows, are advocated to achieve the same effect with significantly less energy consumption in 
driving the flow |^,^ . Much work has been done in this area , but some of the fundamental issues on optimizing 
the design still require a clarification. This paper intends to report some progresses along this line. The emphasis 
on the local Lyapunov exponent of the flow and its role in quantitatively determining the physics of both advective 
'mixing' and diffusive transport, is an unique aspect of this investigation. 

The typical reaction is a three species event, chemical A reacts with chemical B to form chemical C, 

A + B^C. 

The quantities of interest are the concentrations or number densities of each species, Ca,Cb, and Cc- The concen- 
tration of the end product is not independent and can be expressed in terms of the history of Ca and Cb- If the 
chemicals are carried by a flow, the mathematical model for the mixing-reaction process is the advection-diffusion- 
reaction equation. Written out explicitly for each species, they are 

dCA/dt + w^CA=-y ■{DVCa)~kCaCb] (1) 
dCB/dt + V • VCs - V • [DVCb] - kCaCb- (2) 



In the general case, the chemicals react according to rnA + nB IC. The concentrations in equations (|1||^) should 
then be replaced by the scaled concentrations = Ca/itl, Cg = Cb/tl, and = Cc/l- 

We will assume that the original carrier flow is sufficiently energetic that the back-reaction of the reaction process on 
the background flow is negligible. This assumption decouples the third governing equation, for example, the Navier- 
Stokes equation for the flow field, from the two coupled advection-diffusion-rcaction equations (|l|,^. The principles 
for the design optimization, are then found by solving the advection-diffusion-rcaction equations and identifying the 
flow features that crucially affect the production rate and quality. 

The underlying physics is better explained by a transformation of the variables of the original advection-diffusion- 
rcaction equations. Define 4> = Ca — Cb and / = Ca + Cb-, one has 
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9(/)/9t + v- V(/) = V- (DV0), (3) 
df/dt + V • V/ - V • {D\/f) ~ ^K{f - 02). (4) 

Obviously p-(t>'^^ ^CaCb > and / > ||(/)||. 

There are generally two classes of reactors, which we will call a closed flow system and an open flow system. A 
canonical closed flow system is a stirred tank Mathematically it corresponds to an initial value problem, with 
simple boundary conditions. The degree of mixing is measured over time. A canonical open flow system is a tubular 
device such as an automobile catalytic converter. Raw exhaust gas constantly flows in, the poisonous elements are 
removed by chemical reactions, and the treated gas is discharged to the tail pipe. The system is expected to function 
in steady state. Mathematically it is a time-independent boundary value problem. The integral form of equations 
(^,^ makes the global balance transparent. For an engineering device operating in steady state, 

dt dt Jy 

The remaining part is 

J [/V - DVf] -da^-^nj^if- q>^)d'^. 

Here we have assumed incompressibility V • v = for mathematical simplicity. The left hand side is the rate of new 
reactant influx and the right hand side is the rate of rcactant depletion due to reaction. Naturally for a steady state 
device they balance each other out. 

The fast reaction scenario would further reduce the mathematical complexity and help clarifying the basic issues. 
If the chemical reaction rate is sufficiently fast that 

/[/v-^V/]-^a 



then / « \\(f)\\ everywhere in the system since 

/ [f - (/.2]d3x « / fd 
Jv Jv 



3x. 



Mathematically that is to say that the reaction is fast so once (j) is found, the problem is solved by setting ||(/)|| = /. 
Physically it says that two different reactants react so fast that they can not coexist at the same point anywhere in 
the device. Hence at least one of reactant densities has to vanish locally and ||(/>(x)|| = /(x) = Max{CA,CB)- For 
engineering design purposes, the equivalent statement is that the reaction is sufficiently fast so that the efficiency of 
the reactor is determined by the efficiency of mixing. This is precisely the motivation for embedding the reactants 
into a turbulent flow and more recently a chaotic one. We note that the use of a catalyst is intended to boost the 
microscopic reaction rate k. 

As expected, the inequality in equation (^) is equivalent to the statement that the Damkohler number is much 
greater than the Peclet number. The Peclet number Pe is the ratio between the characteristic diffusion time L^/D 
and the advective time scale L/U. The Damkohler number Da is the ratio between the characteristic diffusion time 
and the typical reaction time l/«;(//2) with (• • •) an average over space. A reaction is fast if the typical reaction 
time is much shorter than the advective time scale, i.e. Da/Pe ^ 1. Although in a closed flow system like a stirred 
tank the ratio Da/Pe might be less than one, it has to be greater, and usually much greater, than one in a working 
steady-state device. This is a trivial statement of the fact that the advective time scale L/U in a steady-state device 
is simply the duration for the reactants to stay in the reactor before discharge. It has to be longer than the reaction 
time. Otherwise even a perfectly mixed reactants at the inlet will not form much product at the time of discharge. 
In practice, the requirement for a 'fast' reaction Da/Pe ^ 1 is not difficult to accommodate. For example it can 
always be satisfied in a tubular device by simply increasing the longitudinal length L. Indeed, the minimum design 
constraint of a steady-state chemical reactor is to have the longitudinal Damkohler number greater, or much greater, 
than the longitudinal Peclet number. Here longitudinal refers to the streamwise direction. 

The degree of global mixing in a reactor is measured by 



a(i)^i|^^Vx (6) 
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with 

ip=(j)^ J (f)(fx/ J (fx 

The separation between reactants is always non-negative so cr > 0. The governing equation for the global mixing index 
is 



da 
It 



J'- J D{\7Lpfd^yi. 

The flux of separated chemicals into the system is 

T^- j[^^-DV^]-da.. 

Now we can distinguish two kinds of mixing devices, closed and open, which correspond to an initial value problem 
and a timc-indcpendcnt boundary value problem respectively. For a closed system, 

= 0, da/dt = - y D{Vipfd^x. 

Obviously 

cr > — > (7 = for sufficiently long time. 

An open system is markably different. For example J- normally is positive. It must be positive for steady state 
device where da/dt — 0. Pictorially speaking, separated chemicals are brought in, tf'^ non-zero where v-da < 0. Mixed 
and reacted chemicals are discharged, (/s^ essentially zero where v • da > 0. 

The flow design is driven by the rate of cr ^ which measures the quality of the reactor. The desire for a background 
flow can be clearly illustrated in a closed system, where the relaxation of a is completely due to diffusion driven by 
the gradients of tp. Since 

the relaxation time scale is P/D. If there is no convection, I is of reactor size and the relaxation time scale P/D is 
hopelessly long. 

Towards a better reactor, the carrier flow must be designed to drive gradients in ip or <j), which is an ideal task for 
a chaotic flow. In addition, to achieve better quality of the end product, it is desired that the gradients are driven 
'uniformly' In the case of a chaotic flow, the size of the KAM islands should be minimized, which is an obvious 
direction for optimization. There is another, more subtle, issue concerning the non-uniformity of the chaoticity and 
the existence of a barrier to diffusion even within the chaotic region. This subtlety can not be resolved by the standard 
Poincare section technigue, but it can be made mathematically precise in terms of the finite time or local Lyapunov 
exponents of a flow. The goal of this paper will be to explain how to understand the merit of a reactor by the finite 
time or local Lyapunov exponent of the carrier flow. 

Our approach based on the local Lyapunov exponents should be contrasted with the past heavy reliance on the 
Poincare section as a diagnostic to measure the extent of mixing, an approach that is prone to misinterpretation and 
does not provide any rate information Q. The limitation of the Poincare section technique was explicitly noted by 
Swanson and Ottino in |jll|, and prominently restated by Bryden and Brenner in Muzzio and Liu presented a 
numerical demonstration in a two dimensional mixing-reaction problem The local Lyapunov exponent analysis, 
in our view, complements the Poincare section by providing the exact rate information and removes the ambiguities 
associated with the Poincare plots. 

The rest of the paper is organized as follows. Section gives an introduction on how the finite time Lyapunov 
exponent enters the description of a chaotic flow and the passive scalar transport. Section |l| briefly describes the 
time-dependent solution to the advection-diffusion equation, the more detailed analysis can be found in for the 
two dime nsional c ase and H for the three dimensional case. The main body of the paper, which includes sections 



IV V V] VI] VIII, deals with the practically attractive steady-state reactors. The main mathematical formulation 



has utilized cylindrical geometry and quasi-two-dimensional flows for clarity and convenience, although the most 



important concepts and results can be readily generalized. Section [X highlights the main points of the paper. 
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II. FINITE TIME LYAPUNOV EXPONENT AND A CHAOTIC FLOW 



A chaotic flow is characterized by the exponential variations of neighboring streamlines. For example, the separation 
between neighboring fluid trajectories of a divergence-free, time-periodic two dimensional flow obeys 

{dlf = {dlo ■ e)2 exp(2Ai) + {dk, ■ exp(-2At), (7) 

with A > 0. A flow is called chaotic if 

A°° = lim A > 0, 

t — >oc 

otherwise it is said to be integrable. We note that the laminar to turbulent transition is marked by the emergence 
of many scales for the velocity field. The integrable to non-integrable transition concerns the behavior of the fluid 
trajectories, not the spectrum of the velocity field. This is why the stochasticity of a chaotic flow is also called 
Lagrangian turbulence (Lagrangian trajectories), in analog with the usual Kolmogorov-type Eulerian turbulence 
(Eulerian flow field). As far as the theory of advection-diffusion equation concerns, a qualitative change in the 
solution of the equation occurs at the integrable- nonintegrable transition. This should be contrasted with the standard 
treatment that presumes a turbulent background flow. 

What we mean by a qualitative change in the solution of the advection-diffusion equation can be made precise by 
the rate at which the gradients of passive scalar field grow before the diffusive relaxation dominates the solution. In an 
integrable flow, V(/) grows linearly in time. If the Lagrangian trajectories become chaotic, V0 increases exponentially 
in time. The exponential rate is given by the finite time or local Lyapunov exponent A(^, t), which depends on both 
position and time. That is to say, the value of A depends on where the fiducial trajectory initially starts and how 
long one is tracing the fiducial trajectory. The time and especially the spatial dependence of the finite time Lyapunov 
exponent define the essential characteristics of the passive scalar transport. 

The effects of advection and diffusion can be most clearly seen in Lagrangian coordinates. The Lagrangian coor- 
dinates are defined by a one-to-one mapping between the initial position and the current position of a fiuid element, 
which is found by integrating 

dx{^,t)/dt = v(x, t) 

with the initial condition x(^,i = 0) = ^. The description of a physical phenomena is independent of the choice of 
the coordinate system. Both x and ^ can provide the necessary coordinate system, the first is the usual Eulerian 
coordinates while the second one is the well-known Lagrangian coordinates. Unlike the usual Eulerian coordinates, 
the Lagrangian coordinates have a non-trivial metric tensor. In fact it has both a space and a time dependence. If 
the metric tensor of the Eulerian frame is the unit matrix, the metric tensor of the Lagrangian coordinates is defined 
as 

gij=dx/dC -dx/d^ and 5*^ = Vf • V^^ (8) 

The metric tensor, by its definition, is entirely determined by the fiow field. Even for a perfectly smooth steady fiow 
field, the metric tensor tends to become singular as t becomes large. If the flow is integrable, the metric tensor diverges 
at most quadratically in time |^. For a chaotic fiow, the metric tensor blows up exponentially in some subspaces. As 
long as the fiow field is well-behaved, the metric tensor is well behaved and does not possess any finite time singularity. 

Since the Lagrangian coordinates are attached to the fiuid elements, the solution to the ideal advection equation 
is equivalent to integrating the trajectory x(^,t) and setting (/)(x(^, i), t) — 4>{^,t — 0). If the spatial gradient of 4> is 
written in Lagrangian coordinates, 

[V0(x,i)]' = Vo0(O-?-Vo0(e), 

the coordinate derivative in Lagrangian coordinates of 0, Vo0(C)i would be time independent, i.e. given by initial 
condition, W4>{x,t = 0). When diffusion is included, 4>{£,,t) would have a time dependence, so does Vo4>{£,,t). Two 
important characteristics can be immediately explained once the covariant representation of the metric tensor, g^ , is 
identified as the Oseledec matrix [Q. Since gij is a positive, definite symmetric matrix, it can be diagonalized with 
real eigenvectors and positive eigenvalues, 

gij = exp{2Xt)ee + exp(— 2At)ss (9) 

with A > 0. The gradient of (j) can be written in the form. 
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(V0)^ = (Vo(/) • exp(-2At) + (Vo^ • s)^ exp(2Ai). 



(10) 



The metric tensor that enters the last calculation is its contravariant form g^^ , which is the exact matrix inverse of 
gij, hence g^^ — exp(— 2At)ee + exp(2At)ss. 

What immediately can be seen is that under ideal advection, the passive scalar gradient would blow up exponentially 
in time without bound since Vo4>{0 time-independent. This reflects the singular nature of the diffusion term as 
a perturbation to the ideal advection equation. The existence of a diffusion term, no matter how small D is, would 
eventually remove any gradient in (j). Since the metric tensor has an exponential factor, the diffusion responsible for 
the relaxation of Vo0(^,t) must have a super-exponential dependence in order to overcome the exponential factor in 
the metric tensor. The second note is that the diffusion is essentially one dimensional since the exponential growing 
term in the metric tensor is only in the one dimensional subspace defined by s. 



III. TIME-DEPENDENT SOLUTION TO THE ADVECTION-DIFFUSION EQUATION 

In this section we give an example for a bounded flow, i.e. a closed system. The open flow configuration, the main 
focus of this paper, is treated afterwards. A closed system corresponds to an initial value problem. The equation to 
solve is the advection-diffusion equation, 

dip/dt + vVip = V ■ {DVip). 

We have defined 

if = (j) — (j), and ^ ~ j 4>d'^^- 

Transforming into Lagrangian coordinates, 

{dip/dt)^^\/o-iD*gVo^). 
The tensor diffusivity D = D*g in diagonal form is 

D = De^^^ss + De-'^^^ee. (11) 

We have restricted the discussion to a two dimensional flow. The three dimensional case is treated in [|| . 

The exponential anisotropy of the tensor diffusivity along different directions implies that the diffusive relaxation 
is determined by a one dimensional diffusion equation, 

dip/dt = Vo • De^^'ss • Vq^'- (12) 

The usefulness of above formulation rests on the fact that s converges exponentially to a time asymptotic limit §00 (C)- 
A time- independent vector field §00 (C) allows the usual construction of a coordinate system that makes the calculation 
well-defined. Parameterizing distance along the s line by a scalar field /?, one has 

The exponential convergence rate of s to §00 implies that the calculation based on the time-asymptotic spatial 
coordinates is exponentially accurate. 

The solution to the one dimensional diffusion equation is determined by two quantities. One is the fundamental 
time scale of the problem, the local Lyapunov time of the flow 1/A. The other is the dimensionless number = / 
which is the ratio between the characteristic diffusion time scale /D and the Lyapunov time of the flow. If ^ 1 
which is generally true due to the smallness of the molecular diffusivity, the scalar field undergoes a pure advection 
until time ta = (In2ri)/2A. The ideal advection increases the gradients of the passive scalar field (in Eulerian frame) 
by a factor of f2. There is a rapid diffusive relaxation [super-exponential in dip{^,t)/dp] which removes the gradients 
Viy9(x, t) during a short interval of a few Lyapunov time after ta. 

A clarification should be made on the definition of the Lyapunov time for passive scalar transport studies. To our 
knowledge, in the literature (for a review, see Q), the Lyapunov time of a flow is uniformly associated with the stretch 
rate, i.e. the largest positive Lyapunov exponent of the flow. This is a misconception and the reason can be easily seen 
by examining equations (nPUlCl, |ll|,|l^). The correct statement is that the most negative Lyapunov exponent —A, i.e. 
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the convergence rate, defines the Lyapunov time 1/A for the passive scalar transport. In the case of a two dimensional 
time-dependent divergence-free flow which happens to be the starting point of much of the existing literature, the 
positive Lyapunov exponent has the same magnitude as the negative one due to the constraint of area-preserving, so 
this distinction is not mathematically important despite its physical significance. However, in the more general case 
of a three dimensional flow or compressible flows, a correct understanding of this subtlety is necessary both physically 
and mathematically. 

If the finite time Lyapunov exponent does not vary in space, then one has the ideal situation that the diffusive 
relaxation uniformly removes the gradients. Unfortunately the only flow that is known to have this property is purely 
hyperbolic. An example is the Arnold's cat map. Generic flows are non-hyperbolic and the finite time Lyapunov 
exponent has a peculiar spatial dependence. Roughly speaking, A(^, t) varies smoothly along the s lines, which is 
required for equations (|l2|,^3|) to be well-posed. The variation of A(^, t) across the s lines is pathological: the gradient 
of A(f , t) in directions other than §00 has an exponential dependence in time. The smoothness of A along the s lines 
permits the existence of a class of diffusion barrier in a chaotic region. The pathology of A across the §00 direction 
gives rise to the fractal-like transport in both space and time. 

A crude estimate of the mixing process can be based on the probability distribution function of the finite time 
Lyapunov exponent. If the distribution function is an S function centered at X°° . The relaxation of (V0)^ follows a 
simple route. The gradient grows by a factor of il and is then removed during a Lyapunov time centered on time 
ta = (In2r2)/2A. Equivalently the global mixing index cr(t) stays the same until time ta, after which it exponentially 
decays to zero within a few Lyapunov time. If the finite time Lyapunov exponent has a broadened distribution like 
that shown in Q , the time evolution of the global mixing index can be found by convoluting the ideal mixing curve 
with the A distribution function. For a flow having noticeable amount of remnant integrable region, there is usually 
a significant bump at the left end of the distribution function. Since the distribution function is sampled over the 
entire space, there are actually two types of trajectories which could contribute to this bump. The first is simply the 
integrable trajectory lying on the KAM surfaces. The second one forms the so-called stochastic layer, a sticky region 
surrounding the remnant KAM islands. The bump in the distribution function is responsible for the long tail of the 
a relaxation curve. The experimental observation of a spread in a relaxation is determined by the A profile. 

Although the features of the global solution is roughly given by the mean Lyapunov time, the number associated 
with the local Lyapunov time dictates the local details of the solution. For example, since diffusion occurs only along 
an § line, the place with a peculiarly small finite time Lyapunov exponent would pose as a practical diffusion barrier. 
A more precise description is based on the theory of finite time Lyapunov exponent ]^,p^ . The finite time Lyapunov 
exponent exponentially converges to a form of three parts: 

A(e,t) = ^-f ^ + A-. (14) 

The two convergence functions have entirely different properties. The function A is a smooth function of space and it 
is related to the geometry of the s line by 

§00(0 -VoA + Vo -8,0 = 0. (15) 

The function /(^, i) reflects the pathology in the spatial variation of the finite time Lyapunov exponent. First it does 
not vary along the s direction, 

§oo-Vo/(e,i) = 0. (16) 
Secondly it develops an exponentially growing gradients in the direction across the §00, for example, 

Goo • Vo/ exp(At). 

Of course, the time dependence of /(^, t) has to be bounded by ^/t so limt^oo A — A°°. 

Although equation (na) states that the finite time Lyapunov exponent achieves an extreme when V • Sqo vanishes, 
numerical calculations show that A reaches its local minimum where the s line makes sharp bend. The sharp drop 
in the magnitude of A implies that these sharp bends of an s line are practical diffusion barriers. We note that an s 
line belongs to a chaotic region, and they persist even if the flow is far from integrable. Equivalently speaking, even 
if the flow is driven globally chaotic without discernible KAM islands, there are still regions with retarded mixing. 
Since the production of new chemical C depends on the final diffusive relaxation of A and B, these diffusion barriers 
directly affects the quality of a reactor. 
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IV. STEADY STATE MIXER 



An efficient engineering device employs an open flow operating at steady state. The simplest geometry for such a 
device is a pipe |f7|, ^jT5| so we will consider the cylindrical geometry. The flow field is assumed to be divergence-free 
and quasi-two-dimensional, 

v(a;, y, z) = z X \/il:{x, y, z) -f y)i. (17) 

The absence of a z dependence in Vz{x,y) implies that the flow field is divergence- free in the transverse plane, 
Va:j, • Vxy{x, y , z) = with v^y = z X Wtp(x, y, z). If the mixer is operated at steady state {d4>/dt = 0), the governing 
equation is 

+ ^xy ■ ^xy4> = ^ xy ' D\7 xy(t) + ^D^cf). (18) 

oz oz oz 

We have separated the divergence of the diffusive flux into longitudinal and transversal components. For an operating 
reactor, the longitudinal component becomes exponentially smaller than the transversal component going downstream 
before the onset of diffusive pulse that removes the transversal gradients. In the usual case that D is small and 
B?Vz/LD ^ 1, the longitudinal diffusion term is negligible so 

dcj) 

W^ — -I- Vxy ■ Vxy(j) = ^xy ■ DVxy4>- (19) 

Here R is the radius of the cylinder and L is the total longitudinal length of the device. 
The in-fluxes of reactants A and B are 

Ta,b = J Vz{x,y,z = 0)CA,Bix,y,z = 0)dxdy. 

The flux for each chemical can be defined for any cross section in the reactor, but only their difference is important. 
The differential flux of the reactants 

F(z) = J V;,{x,y)(t){x,y,z)dxdy = J Vz{x,y){CA - CB)dxdy 

is a constant in the longitudinal direction, 

dF{z)/dz = J v.^^^^^^dxdy = 0. 

For practical applications, an initial input is proper if F{z = 0) = 0. The degree of mixing at the transversal plane is 
given by 

a{z) = ^ J Vz{x,y)<j)'^{x,y,z)dxdy. (20) 

This spatially dependent mixing index cr{z) should be distinguished from the global mixing index a{t) defined in 
equation (^ . The reactants become better mixed as they go downstream. The rate of mixing is determined by the 
diffusive effect, 



da 

dz 



J DVxy4> ■ Vxy(j)dxdy. (21) 



The percentage of the chemicals failed to react is roughly measured by 

a(z^ 
aiz = 0) ■ 

The quality of the mixer(reactor) is measured by the Q-factor 

1 ^ a(z^ 
^ 7^ cr(z = L)' 
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In the impractical limit of L — > cx) the chemicals would react completely for a proper initial input F{z = 0) = and 
the Q-factor is infinite. For a finite L, the Q-factor is generally finite as well. For an intake flux 



the amount that does not react is about 



Ta = j Vz{z = 0)CAdxdy, 



If Q is infinite, the output of the final product is 

Oc = J Vz{x,y,z = L)Ccix,y,z = L)dxdy = Ta- 



Otherwise 



The optimization of a mixing device or a reactor is primarily looking for a balance between the production rate 
Ta, reactor length L, and the desired Q-factor. Implicitly through the Q-factor there is also an energy consumption 
penalty (power P) for driving a flow field v(x). These correspond to four design constraints: 1] geometrical constraint, 
i.e. the physical size limit for the device; 2] production rate constraint; 3] production quality constraint; 4] energy 
constraint. The first three are kinematic constraints since they are uniquely determined by the fiow field [Except 
that in a reactor, there is a minimal length constraint due to the finite reaction time, i.e. the system size has to be 
greater than v^/ K{f 12). This point was made earlier in the introduction]. The last one is a dynamical constraint and 
it requires solving the Navier-Stokes equation to relate the driving term to the fiow field. 

This paper addresses the kinematic constraints. The dynamical constraint for the fiow field is more case-dependent. 
Nevertheless, a general statement can be made that it is usually less energy consuming to produce a smooth, non- 
turbulent, but chaotic flow, than a turbulent one. In the case of a tubular device, this is reflected by a small increase 
in the pressure drop over the tube, or equivalently the pumping power, to induce a chaotic flow with satisfactory 
mixing properties [ |l6[ . 

The physical meaning of the kinematic constraints will become clear once the solution to equation ( p^ is found. 
The most important property of the flow field in determining the kinematic constraint is the so-called local Lyapunov 
exponent, just like the time-dependent initial value problem. Instead of time, the z coordinate will be used to 
parameterize a fluid trajectory. The only constraint for z to be a good scalar for parameterization is that nowhere 
vanishes. This is naturally satisfied in our choice of a quasi-two-dimensional flow fleld, dvz/dz = 0. Any point 
satisfying Wz(xo,2/o) = would imply that Vz vanishes along a straight line from {xQ,yQ,z = 0) to (xo,i/0:^ = L), 
which strictly prohibits transport. Hence Vz would be made nonzero in a working device. This is easily achievable by 
an infinitesimal longitudinal perturbation. 



The singular nature of the solution to the equation (19) can be seen by a coordinate transformation. The fiuid 
trajectory is now parameterized by z and given by 

cb ^ v^{x,y,z) ^ dy ^ Vy{x,y,z) ^^^^ 
dz Vz{x,y) ' dz Vz{x,y) 

For an arbitrary initial position {xo,yo) at z = plane, a trajectory of the fiuid element can be traced to(x,y) at 
z-plane. The functional relationship between x(xo,yQ, z) and y{xo,yo, z) is found by integrating equation (G3) with 
initial position (xo,?;o) from z = to z. 

If the mixing equation is transformed into the (a;o,?/o) coordinates, one has 

dcj) 1 

(a^)-- = Vz{xixo,yo,z),yi-o,yo,z))^^^y^ ' ^5V.,,,</>(xo, yo, (23) 

The metric tensor of the (xo^yo) coordinate, 5^, has two forms, covariant and contravariant. The contravariant form 
is defined as 

Va;o • Vyo Vyo • Vyo 

The covariant representation is given by 
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C d{x,y) ^ d{x,y) d{x,y) ^ d{x,y) \ 
dxQ dxQ dxo dyo 1 
dix,y) ^ d{x.,y) d{x,y) ^ d{x.,y) I 
dxQ dyo dyo dyo / 

The metric tensor is a positive, definite matrix, so it can be diagonahzed with positive eigenvalues and real eigen- 
vectors. If the covariant form is written as 

gij — exp(2?7z)ee + exp(— 2?7z)ss (24) 

with r/ > 0, the contravariant form is 

= exp(— 2?7z)ee + exp(277z)ss, (25) 

since gij and g'^^ are matrix inverse of each other. The meaning of e and s can be made precise by following the 
distance between neighboring two points, 

{dxY + {dyf = [dxo,dya]gij[dxo,dyo\^ . 

Hence two initial points separating along the e direction would diverge exponentially going downstream, but they 
would converge exponentially going downstream if their initial separation is along the § direction. Conventionally e 
is called the unstable direction and s is the stable direction. The function ri(xQ,yQ, z) which depends on the initial 
position of the trajectory {xo,yo) and the longitudinal ending point z, is a local Lyapunov exponent measuring the 
exponential separation rate. The inverse of r]{xo-,yo-: z) defines a local Lyapunov scale, i.e. the longitudinal distance 
over which the separation of neighboring trajectories varies by one e-fold. This is also the longitudinal distance over 
which the gradient of a passive scalar field increases by one e-fold in the absence of diffusion. 

The mixing- reaction process is easier to understand if one transforms equation (pl|) into the (a;o,yo) coordinates. 



d<7{z) 
dz 



DVxoya(f'ixo,ya,z) ■ g • Va:oyo(t>{xo,yo,z)Jdxadyo. 



The Jacobian of the (xq, yo) coordinates is identically unity if Vxy is divergence-free. Substituting the diagonal form 
of the metric tensor, equation (Eq), into above equation, one obtains 



da{z) 
dz 



D{[sixo,yQ,z) ■ Vxgyg(j){xQ,yo,z)] exp(2?7z) 



+ [e{xo,yo,z) ■Vx^yg4>{xo,yo,z)] exp{-2rjz)}dxQdyo. (26) 

The reduction in (j{z) can be calculated using a simpler form in the usual case R^v^ti/ D 3> 1. The simplification comes 
from two factors: 1] the diffusive flux is negligible for small z, 2] s(xQ^yQ, z) exponentially loses its z dependence, i.e. 
s(xo,j/Oj^) can be replaced by its z-asymptotic limit §oo(a^O;yo) with an exponentially small correction. Combining 
these, one can calculate the mixing rate with exponential accuracy by 



da{z) 
dz 



-D[s(xo,2/o) ■ Va;oy(,(/)(a;o,yo,^)] exp{2r]z)dxodyo. (27) 



(?)-o,.o = 77Tr^-^^^D^^'^'""'''^^'^^^^i^o,yo,z). (28) 



The gradient of (j){xQ, yo, z) in the {xq, j/q) coordinates is found by solving a simplified form of equation (P3| ) 

1 2r,{xo,yo.z)z_^ 

'dz'-'"'"'' v,{x^,y^,z)dl3 df3 

The P coordinate is defined by a parameterization along the s lines 

dxo/dp cx Soo(a;o, yo) with xo (xq, yo). 

The solution of this one dimensional equation differs from that of the full equation by an exponentially small term in 
z. The discrepancy is negligible if R^Vzij/D ^ 1 and the initial gradient has a scale comparable to R. 

The solution to equations (^,^8|) has remarkable properties for a chaotic flow {t] > 0.) Since the longitudinal 
distance plays the role of time in a steady state system, the usual characteristic diffusion time R^/D gives rise to an 
equivalent characteristic longitudinal length for diffusion R?Vz/D. The ratio of this characteristic longitudinal length 
for diffusion and the Lyapunov length defines 57, 

1]^—^, (29) 
a dimensionless number that is typically much greater than one. 
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V. DEGREE OF MIXING VERSUS LONGITUDINAL LENGTH 



Recall that the transversal mixing index a{z) is defined as 



The variation of a downstream is given by 

da{z) 
dz 



= DV(l){x,y,z) ■V(l){x,y,z)dxdy. 



tj is a monotonically decreasing function which corresponds to the time- irreversibility of the diffusion process. Dividing 
above equation by cr(0) leads to a dimcnsionlcss form 

1 da{z) 2 j D'V(l){x,y, z) ■ V(f){x,y, z)dxdy 



cr(0) dz j Vz(x,y, z — 0)(l)'^{x,y, z ~ Q)dxdy 

In the case that rj is a constant, 



1 daiz) D , e^''^ - 1 , 

exp[ — — + 2r;z]. (30) 



cr(0) dz v^R^ ' 2n 

Equation ( pO| ) can be integrated for an exact analytic expression, 

with Vl defined in equation (|29|). This is a remarkable expression that is best interpreted in the units of Lyapunov 
length. There is little mixing and hence reaction in a region within a critical length 

Lc = \n{2Q.)/2T] (32) 

from the intake of the device. In fact, there is only one e-fold drop in ct(z)/(t(0) over the entire longitudinal length 



ct(O) 



Each additional Lyapunov length beyond produces a super-exponential jump in the production quality, e.g. at a 
length Lc + N/r], 

thus a reator of longitudinal length Lc + N/ -q has a Q-factor of 

Q = exp[e^]. (33) 

The above formula implies that it takes less than = 3 additional Lyapunov lengths to achieve a Q-factor of 10^. 
For a reactor with a length of Lc + 3/77, less than one part of ten thousands input reactants fails to react at the 
moment of discharge. It must be emphasized that is a modest number even if the reactants are extremely difficult 
to mix, which is usually the case due to the smallness of D. For example, even if = 10^°, Lc is about 12 Lyapunov 
length 1/77. A reactor of Q-factor 10^° requires a total longitudinal length of 15 Lyapunov lengths. In giving out these 
numbers, we have assumed that the chemicals undergo a diffusion- limited reaction. Even if the reaction rate is so 
slow that Lc is less than a typical reaction length Vz/ k{J /2)^ an additional longitudinal length of a reaction length in 
the engineering design would provide the extra room to achieve the desired Q-factor. 

The assumption that 77 is a constant, is a deceptively nontrivial one. In actuality it implies a much more stringent 
topological constraint. It is known that only purely hyperbolic system can have a constant Lyapunov exponent, 
the canonical example being the Arnold's cat map. A point is hyperbolic if its stable and unstable directions are 
not degenerate. Generic systems such as hamiltonian flows, are not purely hyperbolic. In fact, we suspect that 
a differomorphism higher that would generally have nonhyperbolic points. In the case of hamiltonian systems, 
remnant KAM tori are one manifestation of the nonhyperbolicity. Points on the chaotic set can also be nonhyperbolic. 
In hamiltonian systems, they are responsible for the sharp bending of the s lines. 
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VI. LOCAL LYAPUNOV EXPONENT AND THE Q-FACTOR 



The local Lyapunov exponent, or equivalently the local Lyapunov length, has a profound role in determining the 
quality of a reactor. In retrospect, the design of a quality reactor can be considered as an optimization of the local 
Lyapunov exponents. The simplest way to understand the effect of a spatially varying Lyapunov exponent on mixing 
is through its probability distribution function P{ri, z ^ L) for a steady state reactor. If 77 is a global constant, 
P{ri, z — L) would, of course, be a delta function. Nonhypcrbolicity prevents this for a generic flow. Instead there is 
a large spread in the distribution function. The mixing index, or the Q-factor, can be approximated by convoluting 
the distribution function with equation (|3^). For a device of longitudinal length L operating with reactants of 
characteristic diffusion time scale /D, there is a critical Lyapunov exponent ?/c(z = L) given by the solution to 
equation 

r ' 1,^!!,^!!^. (34) 



277c (i) D 
As long as 

De 

r]c{L) > rf{L) with [z = L) 



the proportion of the input reactants that fails to react is roughly given by 



assuming that P{7], z — L) la normalized 



Jo 

1*00 

/ P{ri,z^ L)dT] = l. 
Jo 



A crude estimate for the Q-factor of the reactor is then 

rvciL) 



= 1/1 P{v,z = L)dr]. (35) 
^0 



The second critical Lyapunov length 77-^(2 ^ L) — De/2vzP? is determined by 



0. 



Lc is a monotonically decreasing function in 77 if 77 > 77^, but a monotonically increasing function if 77 < rf^ . The 
longer the characteristic diffusion time E? /D, the less likely that the second critical scale rj^{z = L) plays a role. In 
other words, unless E? / D is small enough to render the chaotic mixer unnecessary, one does not need to worry about 
the subtleties associated with rj^ . 

In the language of hamiltonian mechanics, the integrable region or the KAM surfaces, are absolute barriers to 
advective transport because no trajectories could cross them, at least for systems with no Arnold diffusion. It is also 
an effective barrier for diffusive transport in the advection-diffusion theory. This is explained by the time dependence 
of the largest eigenvalue of the metric tensor. It was shown |^ that the largest eigenvalue of the metric tensor in an 
integrable region grows at most quadratically in time, with the prefactor given by the shear rate of the surface 
i.e. the derivative of the rotational transform normal to a KAM surface. Contrasted with the exponential growth in a 
chaotic region, the effective diffusivity in (zq, j/o) coordinates, Dg^^ , is small and hence the diffusive relaxation occurs 
on a much longer time scale. 

The simplest design optimization is to reduce the size of the remnant integrable regions. Once resonant perturbations 
are present, global stochasticity can be achieved by increasing the perturbation strength. This is remarkable since it 
separates the global stochasticity for the trajectories from the spectrum of the velocity field. A turbulent flow field is 
usually associated with a broad spectrum of the velocity field. A laminar, but chaotic, flow has well-behaved velocity 
field, which is reflected by a few isolated spectral peaks. The optimization involves the adjustment of the relative 
strength of the peaks. The reduction of the integrable regions can then be realized under the constraint of maintaining 
the laminar nature of the flow. 

The signature of a significant amount of integrable region in the local Lyapunov exponent is a large bump at the 
near-zero end of the distribution function. The local Lyapunov exponent is also well-defined using the metric tensor 
in an integrable region, even though the asymptotic Lyapunov exponent vanishes. This is trivially explained by the 
fact that a quadratic function can always be approximated by an exponential function locally. 



11 



VII. GEOMETRICAL CONSTRAINT ON A PRACTICAL DIFFUSION BARRIER 



Even if the integrable regions are so small that they are all invisible to naked eyes, or in the extreme, of size smaller 
than DL/vz, the local Lyapunov exponent still has a finite spread in the distribution function. What makes the 
matter worse is that the distribution function tends to preserve an asymmetry biased towards the small 77 end. This 
statement is based on equations ( p"4[p^Jl^ ) and the numerical observation that r] always takes a local minimum, rather 
than a local maximum, where the § line makes a sharp bend. Since the local Lyapunov exponent varies little where 
the s line is straight, but makes a sharp dip when the s makes a sharp bend, there is an overall bias towards small rj 
in the distribution function, even if only the chaotic trajectories are included in the distribution function calculation. 

Surprisingly the presence of a small local Lyapunov exponent does not always lead to a practical diffusion barrier. 
To be a practical diffusion barrier, not only a small effective diffusivity is required, but also the size of the structure. 
Roughly speaking, the diffusivity D sets the minimal size of a practical diffusion barrier, although the integrability of 
the trajectories also plays an important role. In an integrable region, the s lines closes on itself forming closed KAM 
curves. The shear-induced fast diffusion is confined within the KAM surfaces so the radial diffusivity is just D. The 
minimal size requirement of a practical barrier is given by 



w = yjDL/vz- (36) 

It is the mean spread of an initial (5-distribution at the center of an integrable island at the time of discharge. 

In a chaotic region, diffusion occurs only along the § direction, but the s lines can not close on themselves. A single 
s line indeed fills an entire ergodic component densely. The previous estimate for the integrable case, equation (jsF 
usually is too small. The minimal size for the diffusion barrier is comparable to 



i D cxp{2riz)v^^dz. 
Jo 



The exact integration requires an exact form for rj, and the results are much more difficult to interpret. We will use 
an idealized simple form for r] to illustrate some of subtleties. Noting the empirical result that rj cx — (lnK)/z Q with 
K the curvature of the s line bend, we prescribe 

r] = T]c- (co \iik)/z, 

with the understanding that k is scaled by 1/R and hence dimensionless. Combining with equation (^), one finds 

w^R/k"". (37) 

This expression is not valid for large enough k that rj = rjc — (colnK)/z approaching zero. That gives the lowest 
bound w ~ \/dJv^t]c, consistent with the worst scenario predicted by equation (^6|). Otherwise, equation ( ^7|) gives 
a minimal barrier width estimate that depends on the sharpness of the bends. 

The quantity that w should be compared with, is the size of the region having a significantly smaller local Lyapunov 
exponent. This region can be estimated as follows ||^. The position (xq, i/q) with a local minimum for 77 is solved from 

Ke{xQ,yQ) = and W-^ay,, ■ KeixQ^y^) > 0, 

where Ke is the curvature of the e line 

The core of the diffusion barrier is approximately bounded by the curve satisfying 

^xoyo ■ l^ei^Q, I/q) + Kg(a;Q, J/q) = 0. 

(xq, Uq) is the closest point from (xg, y§) satisfying the above constraint. Since the diffusion barrier is a sharp bend of 
an s line, the width of the diffusion barrier is the inverse of the curvature at (xg, j/q), i.e. l/Ke{xQ, i/q). One conclusion 
can be drawn immediately. Those "diffusion barriers" with a size 



l/«e(4>2/o) < VDL/vz (38) 

are not practical diffusion barriers affecting the Q-factor of a reactor. A relaxed criteria is obtained by comparing 
l/Ke(a;Q, ?/§) with w in equation (|37|). This implies that the local geometry of the bends prohibits a simple resolution 
and explicit case-dependent calculation is required. 
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The cut-off size for a practical diffusion barrier brings a natural limit on the finest grid size in a numerical calculation. 
It also opens the possibility for a paradoxical statement that a smaller j] in the distribution function might not imply 
poorer mixing because the geometrical constraint, equations ( |36| , |37| ). Consequently equation ( |35| ) might not be a 
proper estimate for a globally chaotic flow. 

An appropriate formula for estimating the Q-factor is based on the exact spatial dependence of the local Lyapunov 
exponent, ri{xo, yo, L). The formula is based on the integral 

J 6(7;'= - 77(xo, yo, z))dxodyo, 

where Q(x) is the step function 

6(a;) = for a; < 0; 1 for x > 0. 

The geometrical constraint, equation (|3^ ) or alike, is incorporated by imposing a fixed grid size of DL/vz for 
calculating the above integral. Written in discrete form, 

^ = DLj:^j:^U{rf-v{i,j,L)) (^^^ 

where A is the total area of the cross section and the summations are over all grid points on a cross section. In the 
case that 51 3> 1, \/ DL/vz can be tiny compared with the transversal scale of the device R. Equation ( ^9|) is the 
proper formula for calculating the Q-factor in a highly chaotic flow. It should also yield a better estimate for the 
near-integrable cases, but might not be necessary. In principle, if the flow fleld v(a::, y, z) is known, one would calculate 
vixg^ yo, L) and use equation ( ^9| ) for Q. The main usage of a probability distribution function based formula like 
equation (|35| ) is in experimental situations where an exact measurement of rj{xQ,yo, L) is inaccessible. P{ri,z = L) 
can be easily approximated by a limited amount of measurements. For those who could afford the time and resources, 
all quantities can be found to exponential accuracy by solving equations ( |27| , |2^ ) for a given velocity field v{x,y, z). 

In retrospect, a fully chaotic reactor can be further optimized by balancing the value of the local Lyapunov exponent 
and the sharpness of the s-bend. The basic idea is to either straight out the s line or make the bends sufficiently 
sharp that equation (|3^) is violated. 



VIII. A SPECIFIC NUMERICAL EXAMPLE 



Since the goal of this paper is to demonstrate the basic principles rather than working with a specific device, the 
illustrative example will be chosen as simple as possible and computationally as efficient as possible. Nevertheless, 
there are a few features that we do intend to include to make it practically relevant. First we want the example 
flow to have a bounded chaotic region since it is supposed to be confined by a pipe. The flow is also expected to be 
stochastic along the axis. Divergence-free is another useful feature to be consistent with earlier analysis. 

The simplest chaotic flow of the form equation has Vz a constant and tp a periodic function of z. The standard 
treatment of this class of flow is to construct a mapping by sampling at the period of ip along z axis. The mapping 
thus generated preserves both the topology of the Lagrangian trajectories, the functional form of the local Lyapunov 
exponents (except z takes discrete values), and the geometry of the s lines. Mapping is enormously more efficient 
than the flow in computations. Hence our example will employ a mapping directly, but with the understanding that 
it was reduced from a constant-?;^ and z-periodic-'0 smooth three dimensional flow. 

It should be noted right away that in a practical application, an analytic form of the flow field is unlikely available, 
let alone an exact reduction to a mapping. However, this simplification is just for the convenience of illustration, the 
lack of it does not prevent a practical calculation. 

There are three equations mapping ?/„, z„) to {xn+i, yn+i, -Zn+i)- One of them is simply 

Zn+l = Zn + Az 

with Az the period of ip. To insure area conservation, the example mapping that relates {xn,yn) to {xn+i,yn+i) is 
defined by a generating function, which has one free parameter fc, 

S{Xn,yn+l) = Xnyn+1 + fcln(l + X^ + y^+i). 

Here k signifies the perturbation. The map is given by 
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FIG. 1. Poincare plot of the map given in equations with k — 1.4. The solid line is an s line on {xo,yo) plane. 



2fc , 7V 2 ; (40) 



dS _ x_ 

dxji 1 -|- x^^ -\- y^-^j^^ 



Vn^ — ^ Vn+l + , o " „2 ■ (41) 



For k <ti 1 this is the twist map with the angle of twist per iteration or every longitudinal advance of Az, equal to 
2fc/(l + a;^ +2/^)- The map has a non-zero Lyapunov exponent near the axis, ^ 1, if is greater than unity. 

This map is made explicit by solving a cubic equation for in terms of a;„ and y„. To do this, define s so 

Un+i = s + j/„/3, which implies y^^i = + 2sy„/3 + j/^/9. One then finds 



which implies 



Define 



2 2 J^oo 



+ s(l + xl + _ 4^y2 ^ 2/ca;„ - ^y„(l + xl + iyfj = q. 



1 2 1 

a = 1 + 1 + a;^ - -y,^ and b = 2fca;„ - -?/„(l + a;^ + -y^^J. 



The equation for s can then be written in the standard form for a cubic, s'^ + as + b — 0. This equation has an unique 
real root if 

is positive, the situation for this map. The real root of the cubic is then 

s = (r- -)3 - (r + -6)3. 

It is straightforward to show that the map is stochastic near the axis for fc > 1. Near the axis, x^^ + << 1, the 
map reduces to yn+i = y-n — 2fcx„ and Xn+i = 2;„ + 2fc?/„+i. This map is linear and can be solved by Xn = xoF" and 
y„ = yoF". One finds (F - l)j/o = ~2kxo and (F - 1)xq = 2kTyo or (F - 1)^ + {2k)^T = 0. So 



F = (1 - 2fc2)+2v/fc2(fc2 - 1). 
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Since the map is area preserving the two roots must satisfy r+r_ = 1. If fc > 1, one of the roots, r_ , satisfies |lr|l > 1 
and the map is stochastic with a Lyapunov exponent In ||r_||. 

Figure 1^ plots the intercepts of trajectories on the {xq, yo) plane. The topology of this intercepts tells the integrability 
of the trajectories going downstream. Since the flow field is periodic in z, the cross sections at z — NAz would have the 
same topology as that of (a;o, yo) plane. Because of this, an efficient way to understand the topology is by projecting 
the intercepts of a long trajectory with cross sections at z = NAz,N — 0, 1,2,---, on the {xQ,yo) plane. Since 
we already have the proper mapping, this becomes simply plotting the trajectories of the two dimensional mapping 
y„+i). This is exactly how figure is generated. Obviously our example map has a bounded chaotic 
region and it is stochastic near the axis. 

The most important feature we wish to demonstrate here is the geometry of the s lines and the spatial variation 
of the local Lyapunov exponent. The s(a;o,?/o) line sits on the z = cross section. The statement that diffusive 
relaxation only occurs along the s lines in {xo,yo) coordinates is equivalent to say that, initial chemical reactants 
located along an s line will come close together and react to form new product Cc- The sharp bends of the s line 
give rise to a peculiarly small local Lyapunov exponent which hinders diffusive relaxation. This relationship between 
the local Lyapunov exponent and the s-bends is shown in Figure ^. Wherever the curvature of the § goes up sharply, 
there is a significant dip in the local Lyapunov exponent. 



IX. CONCLUSIONS 



We have considered the design criteria of a chemical reactor device based on a chaotic flow. Particular attention 
has been paid to a steady state reactor, where the Lyapunov length provides the fundamental spatial scale for the 
problem. The advantages of a chaotic flow for a steady-state reactor can be summarized into two scaling relations: 1] 
the minimum longitudinal length of a reactor, Lc, has a logarithmic dependence on the diffusivity D, equation (|3^); 
2] the quality of the reactor, the so-called Q-factor, has a super-exponential dependence on each additional Lyapunov 
length beyond Lc, equation (33). 

The Lyapunov length l/rj, as defined through equations ( p^ , p4[ ), is a local Lyapunov length since it depends on 
the initial position (xq, yo) and the longitudinal displacement z of the test fluid point. The local Lyapunov length are 
defined everywhere in the (xq, yo) plane, independent of whether the trajectory is chaotic or integrable. If the global 
flow field is known, the Q-factor of the reactor can be calculated via equation (39). In the experimental situation 
that the global flow field is not available, one can approximate the probability distribution function of the local 
Lyapunov exponent (length) by repeated measurements over time. The Q-factor can then be estimated by equation 
( ^ ) . Although the derivation is presented with the assumption that flow field obeys equation ( pT| ) , these formulae are 
expected to be useful under more general conditions, such as a weak z-dependence of Vz and a weak time dependence 
of V. The fast reaction scenario, equation (^), which brought mathematical simplicity to the discussion, can also be 
relaxed. A naive way to approach this is assuming that there is no reaction till length Lc + N/r], after which the 
reaction is suddenly turned on. If the fast reaction condition is satisfied, near-perfectly mixed chemicals A and B 
immediately react to form C. Otherwise there is a delay to achieve the same Q-factor. The additional length for the 
reactor is the reaction length, the product of the reaction time 1/k(//2) and the longitudinal flow velocity v^. 

For a closed flow reactor such as a stirred tank, the results obtained for the steady-state device can be straight- 
forwardly translated by identifying the longitudinal length z as time t. The time scales replace the spatial scales as 
the quantities of concern. For example, the Lyapunov length l/jy is replaced by the Lyapunov time 1/A and the 
crtical length Lc is replaced by a critical advection time ta- One difference is that in a steady-state tubular device, 
the streamwise distance can always be used to reduce the problem to a two dimensional one, equation (p2[), but the 
time-dependent closed flow reactor usually has to deal with a truely three dimensional flow, in addition to a possible 
time dependence. Fortunately, aside from some additional subtleties, the main physical charatersitics of the transport 
of a passive scalar in a three dimensional flow is the same as the two dimensional case H . 
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FIG. 2. The local Lyapunov exponent ai z = NAz with N — 5 {Az is the spatial period of the flow field in z direction) and 
the curvature k of the s line are plotted along an s line shown in figure 0. The local Lyapunov exponent makes a sharp dip 
wherever the § line makes a sharp turn. The missing dips at extremely large curvature are due to finite machine precision and 
integration stepsizes. 
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